1. Field of the Invention
The present invention relates to methods of obtaining uniform and independent random numbers on computers. More specifically, said methods relate to the multiplicative congruential generator comprising                1. a positive integer d called modulus,        2. a positive integer z (0<z<d) called multiplier that is coprime with d, namely that shares no common prime factor with d, and        3. a positive integer n (0<n<d) called initial value or seed coprime with d, generating a sequence of integers {r1, r2, r3, . . . } by recursive congruence relationsr1=n, rj+1≡zrj mod(d), 0<rj<d, j=1, 2, 3, . . . ,and giving an output sequence {v1, v2, v3, . . . } of rational numbers in the interval (0,1) byvj:=rj/d, j=1,2,3, . . . ,In contrast to the prior art that takes the modulus d for an odd prime p or a power 2i for a large index i, the invention selects those d comprising two odd distinct primes as the most appropriate, based on the survey of conceivable candidates of composite d. With other key specifications for the choice of said odd primes and for due restrictions on the design of the multiplier z, the invention will render forces to equip any computers and any operating systems feasibly with generators of uniform and independent random numbers with long periods and tested statistical high quality without improbable symmetry restrictions.        
2. Description of the Related Art
Random number sequences on computers are required to have reproducibility and transportability for debugging simulation programs. Said reproducibility is the nature that the identical sequence of random numbers may be generated by demands of users. Said transportability refers to the generation of one and the same sequence of random numbers on transplantation to any other computers or computer languages.
Since the vast amount of random numbers needed for simulations never allows their storage in memory, the sole possible means left is some computational procedure to generate random numbers sequentially. Said requirements of reproducibility and transportability restrict the procedure to the arithmetic of integers which by nature is devoid of truncation and round-off errors on any computers and in any computer languages. Hence a uniform random number sequence on computers should exclusively be realized by some arithmetic as a sequence {x1, x2, . . . } of integers in a range [0, z) for a large but finite integer z, and should be transformed as {u1=x1/z, u2=x2/z, . . . } to the output sequence of rational numbers in the interval [0, 1) by division at every use.
For arbitrary positive integers z and T, any sequence {x1, x2, . . . , xT} of zero or positive integers in [0, z) may be concatenated to form a base-z numeral sequence,x=0·x1x2 . . . xTx1x2 . . . xT . . . =x1/z+x2/z2+x3/z3+ . . .with period T. Being excluded of uninteresting cases that said integers {x1, x2, . . . , xT} are all zero or all z:=z−1, said x is a rational number in (0, 1) and has the expressionx=(x1zT-1+x2zT-2+ . . . +xT)/(zT−1)=n/d, where and hereafter n/d stands for an irreducible fraction with 0<n<d. Since d divides zT−1, d and z are coprime. Any length-T sequence {x1, x2, . . . , xT} of integers in [0, z), which may be an output of any random number generator on a computer, or a length-T sequence taken from natural or physical uniform random numbers discretized with the precision 1/z, thus has an obvious correspondence to a period of the numeral sequence to base z, representing an irreducible fraction n/d with z and n coprime to d.
The division process of the irreducible fraction n/d with 0<n<d to the base z is expressed in the identity and the inequalityrj+1=zrj−xjd, 0<rj<d, j=1, 2, 3, . . . ,where {r1=n, r2, r3, . . . } are remainders of division that never vanish since n and z are both coprime with d. Upon division by zd these give a crucial estimate found by Nakazawa (2004),0<rj/d−xj/z=rj/d−uj=(rj+1/d)(1/z)<1/z, uj=xj/z, j=1, 2, 3, . . . .Thus, numbers {u1=x1/z, u2=x2/z, u3=x3/z, . . . } in [0, 1), irrespective of whether they are created by some computational methods or obtained from any physical random numbers, are approximated from above within a uniform error bound 1/z, by numbers of the sequence {v1=r1/d, v2=r2/d, v3=r3/d, . . . } formed by remainders of the base-z division process of an irreducible fraction n/d, if the denominator d is allowed to be any positive integer to which the base z is coprime.                Nakazawa (2004)/H. Nakazawa: Coset representations of uniform and independent random number sequences I—Cosets associated with periodic sequences. (URL) http://www10.plala.or.jp/h-nkzw (04eprsq1.pdf/Jul. 16-Aug. 12, 2004).        
The term rj in the sequence {r1=n, r2, r3, . . . } of remainders of base-z division n/d is congruent to nzj−1 modulo d for j=1, 2, . . . , with the value of rj taken in (0, d). If the modulus d is an odd prime or a large power 2i of 2, the sequence {v1=r1/d, v2=r2/d, . . . } is the multiplicative congruential random numbers of the prior art generated by the multiplier z, for the initial value (or seed) n and in the modulus d. We analyze the most general form of the modulus d characterized by its factorization into prime factorsd=p1i1p2i2 . . . psis, ij≧1, 1≧j≧s, in order to find and exploit the optimal technological setting.
Let d≧2 be any positive integer, and denote Z*d for the set of integers coprime to d. Equivalence or identity of elements of Z*d is defined by the congruence modulo d; constituents of Z*d may thus be represented by integers in the interval (0, d), which will be replaced with [0, d) for later conveniences. The set Z*d forms a group with respect to the multiplication modulo d, and is called the reduced residue class group modulo d. The order, or the number of distinct elements, of the group Z*d is defined to be Euler's function and denoted as φ(d). For a prime d=p, Z*p consists of {1, 2, . . . , p−1}, giving φ(p)=p−1. For a power pi (i≧2) of a prime p the group Z*pi may be taken to consist of integers in the interval [0, pi) with the exclusion of multiples of p,0=0·p,1·p,2·p,3·p, . . . ,(pi-1−1)·p, totaling to pi-1 in number. These give the following fundamental formula for any prime p:φ(pi)=pi−pi-1=pi(1−1/p), i≧1,which is always even with the sole exception of φ(21)=1. For any z∈Z*d define the cyclic sequence <z>, i.e. the cyclic (sub)group generated by z denoted in sequential form<z>:≡{1=z0,z1,z2, . . . } mod(d).With all constituents taken in the interval [0, d), <z> is identical to the multiplicative congruential sequence generated by the multiplier z with the modulus d and the initial value n=1. For the initial value n∈Z*d the multiplicative congruential sequence is the coset sequencen<z>≡{n,nz,nz2, . . . }≡{r1,r2,r3, . . . } mod(d),again with all constituents taken in the interval (0, d). Periods of <z> and n<z> are one and the same, to be denoted T. Lagrange's theorem proves that T is a divisor of the order φ(d) of Z*d. If the group Z*d is itself cyclic and z is its generator, then Z*d is exhausted by the sequence <z> with the full period T=φ(d). Any coset n<z> is then identical with <z> as a periodic sequence except for the starting element n. The group Z*d is known to be cyclic only for following cases of the modulus d:                1. d=2, 4; cases d=2i, i≧3 are excluded.        2. d=pi with an odd prime p and any integer i=1, 2, . . . .        3. d=2pi with an odd prime p and any integer i=1, 2, . . . .Case 1 is not interesting and discarded. Case 3 will find its natural treatments in the decomposition of the modulus d into prime factors.        
The group Z*d with composite modulus d formed by coprime factors acquires a significant structure regulated by Chinese remainder theorem. Fix an arbitrary integer s≧2 and take pairwise coprime integers d1, d2, . . . , ds. Define d:=d1d2 . . . ds, andej:=d/dj, 1≦j≦s. Since ej consists of dk other than dj, dj and ej are coprime. Euclidean algorithm gives integers Dj and Ej that yield the identityGCD(dj,ej)=1=djDj+ejEj, 1≦j≦s. Define Uj:=ejEj determined by d1, d2, . . . , ds. Since ej consists of dk's other than dj, the above identity gives, with δjk for Kronecker's delta,Uj≡δjk mod(dk), 1≦j, k≦s. Here and hereafter 1≦j, k≦s will be used for short of 1≦j≦k and 1≦k≦s. Take arbitrary integers {zj|1≦j≦s}. The set of congruence relationsz≡zj mod(dj), 1≦j≦s, have a solution z=z1U1+ . . . +zsUs. If z′ is another solution of the above congruence relations, then z−z′≡0 mod(dj) holds for 1≦j≦s, implying z−z′ is a common multiple of d1, d2, . . . , ds. Hence z−z′ is divisible by their least common multiple d, proving that the solution of simultaneous congruence relations is unique modulo d with the formz≡z1U1+z2U2+ . . . +zsUs mod(d=d1d2 . . . ds)for integers U1, U2, . . . , Us determined by d1, d2, . . . , ds. Since this is a congruence relation modulo d, integers U1, U2, . . . , Us may be replaced with any set of integers equivalent to them modulo d. This is a form of the Chinese remainder theorem that holds true for any number s=2, 3, . . . of coprime factors.
Let again d1, d2, . . . , ds be pairwise coprime positive integers for any s≧2, and set d=d1d2 . . . ds. Any element z∈Z*d is coprime to any of d1, d2, . . . , ds, so that there is an element zj∈Z*dj giving z≡zj mod(dj) for each j, 1≦j≦s. Chinese remainder theorem thus proves that this element z taken arbitrarily in Z*d has the expressionz≡z1U1+z2U2+ . . . +zsUs mod(d=d1d2 . . . ds),or that any z∈Z*d corresponds to a unique set (z1, z2, . . . , zs) which is an element of the product set Z*d1×Z*d2× . . . ×Z*ds. Conversely, any set of integers(z1,z2, . . . ,zs)∈Z*d1×Z*d2× . . . ×Z*d,corresponds to simultaneous congruence relationsz≡zk∈Z*dk mod(dk), 1≦k≦s, which determine an element z∈Z*11d2 . . . ds in the form z≡z1U1+z2U2+ . . . +zsUs mod(d) uniquely modulo d=d1d2 . . . ds by Chinese remainder theorem. These establish the one to one correspondence between all elements of Z*d1d2. . . ds and all elements of the product set Z*d1×Z*d2× . . . ×Z*ds. Counting elements in two sets, the multiplicative formula of Euler's function for pairwise coprime integers d1, d2, . . . , ds is obtained:φ(d1d2 . . . ds)=φ(d1)φ(d2) . . . φ(ds).
As already noted, a composite modulus d for the group Z*d must not include any power 2i (i≧2) of 2 in order for <z> and n<z> in Z*d to realize uncorrelated sequences of random numbers. See Nakazawa and Nakazawa (2008a,b) for this discovery; Detailed Description of the Invention also gives related inferences with FIGS. 7A and 7B. Therefore, the most general form of a modulus d relevant to uniform and independent random number problems is summarized in the decompositon to prime factors asd=d0d1d2 . . . ds, s≧1, dk=pkik, ik≧1, 1≦k≦s. Here p1, p2, . . . , ps are distinct odd primes, and d0 is either 20=1 or 21=2. The reduced residue class group Z*2={1} is a trivial group and φ(2)=1. The multiplicative formula of Euler's function gives the order of the group Z*d for the above noted modulus d asφ(d)=φ(p1i1p2i2 . . . psis)={p1i1-1(p1−1)}{p2i2-1(p2−1)} . . . {psis-1(ps−1)},for all cases of our interest.                Nakazawa and Nakazawa (2008a)/H. Nakazawa and N. Nakazawa: Designs of uniform and independent random numbers with long period and high precision—Control of the sequential geometry through product group structures and lattice configurations—. (URL) http://www10.plala.orjp/h-nkzw (03090403.pdf/Mar. 9-Apr. 2, 2008).        Nakazawa and Nakazawa (2008b)/H. Nakazawa and N. Nakazawa: Designs of uniform and independent random numbers with long period and high precision—Control of the sequential geometry through product group structures and lattice configurations—(URL) http://www10.plala.or.jp/h-nkzw (3978erv.pdf/Mar. 9-Jul. 8, 2008).        
Two or more odd prime factors in the composite modulus d introduce into the group Z*d new structures that are absent in moduluses formed by a single prime or by its power. Take for s≧2 the modulus d=d1d2 . . . ds with pairwise coprime positive integers d1, d2, . . . , ds, each of which is assumed to give a non-trivial group Z*dk for 1≦k≦s. Any elements z, n∈Z*d determine corresponding elements in Z*dk byz≡zk mod(dk), n≡nk mod(dk), 1≦k≦s. These give simultaneous congruence relations nzj−1≡nkzkj−1 mod(dk) for 1≦k≦s, which have the solution by Chinese remainder theorem in the following form:nzj−1≡n1z1j−1U1+ . . . +nszsj−1Us mod(d=d1d2 . . . ds), j=1,2,3, . . . ,As noted, integers U1, U2, . . . , Us are determined by {dk|1≦k≦s} and independent of {nk, zk|1≦k≦s}. The sequences <z> or its coset n<z> in Z*d1d2 . . . ds is thus a shuffling, in the sense of random number problems, of component cyclic sequences <zk> or cosets nk<zk> for 1≦k≦s, respectively. Note that this shuffling, though misstated to be linear in Nakazawa and Nakazawa (2008a) based on the appearance of the r.h.s. of nzj noted above, is nonlinear in its essence. This is a congruence relation and the value of the l.h.s. jumps nonlinearly as the r.h.s. crosses any multiple of d=d1d2 . . . ds. See Nakazawa and Nakazawa (2008b). Features of the generator with composite modulus d=d1d2 . . . ds thus stem from those in component groups {Z*dk|1≦k≦s}. Take typically the period of the coset sequence n<z>, including the case n=1 of cyclic sequence <z>. The expression of nzj−1 noted above, which is unique modulo d, proves that the period T of any coset n<z> is realized if and only if all component cosets {nk<zk>|1≦k≦s} resume their respective initial values simultaneously. Let nk<zk> have the period Tk for 1≦k≦s. The period T (or the order) of <z> and n<z> is thus given by the least common multiple,T=LCM(T1,T2, . . . ,Ts).
Resume the most general composite modulus d of our concern with the formd=d0d1d2 . . . ds, d0=1 or 2, dk=pkik, ik≧1, 1≦k≦s, s≧2where p1, p2, . . . , ps are distinct odd primes. We now specify the largest period attainable with this modulus. The trivial group Z*2={1} for d0=21 always has the period T0=1 and, together with the case d0=20=1, Z*d0 is irrelevant in problems of the period. Lagrange's theorem states that the period Tk of nk<zk> for any elements zk, nk∈Z*pkik is a divisor of(the order of the group Z*pkik)=φ(pkik)=2pkik-1qk, qk:=(pk−1)/2, 1<k≦s; note that pk is odd and pk−1 is even. Therefore, the period T of the coset sequence n<z> in the group Z*d satisfiesT=LCM(T0,T1,T2, . . . , Ts)≦LCM(1,2q1p1i1-1,2q2p2i2-1, . . . ,2qspsis-1).Since the group Z*pkik is cyclic with generating elements for 1≦k≦s, the summary of related facts given below in Theorem A is now obvious                Theorem A Let the modulus d be given with s≧2 and with distinct odd primes p1, p2, . . . , ps as d:=d0d, d:=d1d2 . . . ds, dk=pkik, ik≧1, 1≦k≦s, where d0 is either 20 or 21 and gives φ( d)=φ(d). Define integersqk:=(pk−1)/2≧1, 1≦k≦s, with φ(dk)=φ(pkik)=2qkpkik-1 for 1≦k≦s. For any elements z, n in the group Z* d the period T, which is common to the cyclic sequence <z> and its coset sequence n<z>, is a divisor of φ(d), and fulfills inequalitiesT≦2p1i1-1p2i2-1 . . . psis-1LCM(q1,q2, . . . ,qs)≦φ(d)/2s-1.The largest possible period is φ(d)/2s-1. It is realized when d0=20=1 and the following set of conditions (I) and (II0) are satisfied, or when d0=21=2 and the set of conditions (I) and (II1) are satisfied:        (I) Distinct odd primes p1, p2, . . . , ps are chosen so as to give pairwise coprime integers q1, q2, . . . , qs.        (II0) For d0=20=1 with d=d the multiplier z is defined by generating elements zk of the component group Z*dk=Z*pkik for 1≦k≦s asz≡zk mod(dk), 1≦k≦s         in the form z≡z1U1+z2U2+ . . . +zsUs mod(d).        (II1) For d0=21=2 with d=2d the multiplier z is likewise defined by generating elements zk of the component group Z*dk=Z*pkik for 1≦k≦s asz≡1 mod(2), z≡zk mod(dk), 1≦k≦s         in the form z≡1·d+(d+1)(z1U1+z2U2+ . . . +zsUs)mod( d=2d).For the multiplier z realizing the largest period, the group Z* d is divided precisely into 2s-1 portions by the cyclic sequence <z> and its cosets.        
The shuffling arising with moduluses formed by two or more odd prime factors induces effects reflecting structures of component groups, as is the case with periods noted in Theorem A. However, some drastic changes may also arise with the shuffling. The plot of the present invention is to make essential use of such a mechanism, which has its root in an obvious property of cyclic sequences and their cosets in reduced residue class groups.                Lemma B Let d>4 be any integer, and z and n in (0, d) be any elements of the group Z*d which is not assumed to be cyclic. Denote the coset sequence n<z> as,n<z>={r1,r2, . . . }, <rj<d, r1=n, rj≡nzj−1 mod(d), j=1,2, . . . ,If the cyclic sequence <z> includes −1≡d−1 mod(d) as an element, the order T of n<z>, the number of elements or the period, is even, and any halves of the coset sequence are completely correlated as follows with q=T/2,C+:={r1,r2, . . . ,rq}≡{n,nz, . . . ,nzq−1} mod(d),C−:={rq+1,rq+2, . . . ,r2q}≡{nzq,nzq+1, . . . ,nz2q−1}≡{−n,−nz, . . . ,−nzq−1}≡−C+mod(d).        Corollary to Lemma B Let d>4 give a cyclic group Z*d with a generator z. Then an arbitrary integer n in Z*d gives a coset sequence n<z> sweeping all elements of Z*d, but its first and second halves C± are completely correlated. Therefore, n<z> may be used, if at all, only for the half period T/2=φ(d)/2 for a sequence of independent random numbers.(Proof) Let z have the order or the period T giving zT≡1 mod(d). By the assumption of −1∈<z> there is an integer q in 1<q<T≦φ(d) that gives −1≡zq, z2q≡1 mod(d). Thus, 2q=T, T is even, and the rest of Lemma B is obvious. By −1≡d−1 mod(d) the group Z*d invariably includes −1. Though a cyclic sequence <z′> in Z*d may be a proper subgroup and −1 need not necessarily be included in it, the assumption of Corollary gives −1∈Z*d=<z>, proving that Lemma B holds true with <z>. ▪Primitive root generators have frequently been warned that only small portions of periods should be used for random numbers. Irrespective of whether a group Z*d is cyclic or not, its arbitrary cyclic sequence <z> with −1=zq mod(d) should be warned similarly. We may proceed a step further to the original sequences of integers corresponding to cosets C±, namely sequences corresponding to irreducible fractions n/d and (d−n)/d=1−n/d, respectively. Numeral expressions of these fractions to base z,n/d=0·x1x2 . . . xqxq+1xq+2 . . . x2q . . . ,1−n/d=0·xq+1xq+2 . . . x2qx1x2 . . . xq . . . ,reveal simple relations xj+xq+j=z−1= z for all j=1, 2, . . . , and show how the correlation between C+ and C− of the cyclic sequence <z> including −1 reflects the original structure in the sequence {x1, x2, . . . } of integers with improbable symmetry between its first and second halves.        
Take for a while the simplest composite modulus d=p1p2 formed by two odd distinct primes p1 and p2. The largest order of a multiplier z, or the largest period of <z> on this modulus d=p1p2 was shown to be T=φ(d)/22−1=φ(d)/2, which is realized if primes p1, p2 are so chosen as to give coprime integersq1:=(p1−1)/2, q2:=(p2−1)/2;and if the multiplier z∈Z*p1p2 is constructed with primitive roots z1, z2 of moduluses p1, p2, respectively. Note that the coset sequence n1<z1> for any n1∈Z*p1 has the excellence that the sequence sweeps almost all numbers in [0, p1), realizing a nearly perfect uniformity in one period. The matter is the same with n2<z2>. However, there is a spell of Corollary to Lemma B that each of n1<z1> or n2<z2> in its solitary use should be restricted to halves of their respective periods, thus spoiling their noted excellences in uniformity. It is to be admired that the shuffling that works in the composite modulus d=p1p2,nzj−1≡n1z1j−1U1+n2z2j−1U2 mod(d=p1p2), j=1, 2, . . . ,breaks this spell on component sequences. First, they are allowed to repeat their periods any number of times, and problems of (repeated appearances of) their correlated first and second halves are dismissed. Second, the shuffling gives chances to construct the coset n<z> that is free from the similar spell. The following Theorem C, which is the heart of the present invention, clarifies this magical work of Chinese remainder theorem. This is a generalization of the result obtained in Nakazawa and Nakazawa (2008a,b).                Theorem C Let the modulus d be given with s≧2 and with powers of distinct odd primes p1,{grave over (p)}2, . . . , ps as d=d0d, d=d1d2 . . . ds, dk=pkik, ik≧1, 1≦k≦s, where d0 is either 20 or 21. As in Theorem A define integersqk:=(pk−1)/2, 1≦k≦s, assuming that primes p1, p2, . . . , ps are so given that q1, q2, . . . , qs are pairwise coprime. Take likewise a generator zk in the component group Z*dk for 1≦k≦s and construct a multiplier z∈Z*d with congruence relations z≡zk mod(dk) for 1≦k≦s to obtainzj≡z1jU1+z2jU2+ . . . +zsjUs mod(d), j=1,2, . . . ,realizing the largest period T=φ(d)/2s-1 for <z>.            (I) If integers q1, q2, . . . , qs consist of one even and all other odd integers, the cyclic subgroup <z> does not include −1.    (II) If the remainder of the dichotomy holds true, i.e. if q1, q2, . . . , qs are all odd, then the cyclic sequence <z> includes −1.(Proof) Assume the case d0=1 with d=d. Take a component cyclic group Z*dk (1≦k≦s) with dk=pkik and its arbitrary generator zk. This group contains −1≡dk−1 mod(dk) and has the order φ(dk)=pkik−1(pk−1)=2qkpkik−1. Therefore, zkj≡−1 mod(dk) arises when and only when j is an odd multiple of tk:=φ(dk)/2=pkik−1 qk. Since zj≡1 mod(d) occurs for the first time in j>0 at j=T=φ(d)/2s-1, we have zq≡−1 mod(d) at q=T/2 withq=T/2=φ(d)/2s={φ(d1)/2}{φ(d2)/2} . . . {φ(ds)/2}=t1t2 . . . ts.There is also an expressionzq≡z1qU1+z2qU2+ . . . +zsqUs mod(d),which means that zq≡−1 mod(d) can arise only when zkq≡−1 mod(dk) holds for all 1≦k≦s. Therefore, −1∈<z> is true if all of factorstk=φ(dk)/2=pkik-1 qk, 1≦k≦s, are odd, and false if one or more of them are even. Since all of primes p1, p2, . . . , ps are odd, the former corresponds to the case (II). The latter can arises only as the case (I), because q1, q2, . . . qs are restricted to be pairwise coprime. As regards the case d0=2 with d=2d, the largest period multiplier takes the formzj≡1·d+(d+1){z1jU1+z2jU2+ . . . +zsjUs} mod( d=2d);here zk, Uk for 1≦k≦s are those given in the preceding case d0=1. Note that d+1 is even, and zj≡1≡−1 mod(2) holds for any j. Thus zj≡−1 mod( d=2d) realizes if and only if some j gives zkj≡−1 mod(dk) for every 1≦k≦s, which was proved above to be the case in (II), and not the case in (I). ▪The conclusion of Theorem C reveals that an appropriate choice of the composite modulus d, namely            among pairwise coprime integers q1, q2, . . . , qs one is even and all others are odd,breaks the spell that a largest period multiplier z is forced to generate a sequence <z> whose every pair of halves are stringently tied by symmetry. It will be notable that this is accomplished magically by using component sequences <zk> under the same spell constructed on a generator zk of the group Z*pkik. Further implications of this significant comprehension will be seen below in intuitive, geometrical terms.        
Spectral tests have been the most powerful tool to assess the performance of a multiplier z of multiplicative congruential generators, typically with the modulus d formed by an odd prime d=p, or by a power of 2, d=2i. They remain to be an indispensable and powerful weapon to test multipliers of composite modulus generators, hence to test any random number sequences of finite length for their uniformity and independence, if they are computable at all. The inclusion of composite modulus cases, however, necessitates many reflections on structures to grasp correctly what has so far been passed unnoticed. The modulus in this paragraph is any positive integer d>4. The multiplier z and the initial value n are any elements in the group Z*d. The group Z*d need not be cyclic, nor z need be a generator. Let l=2, 3, . . . be an arbitrary integer, and let the coset sequence n<z> have the period T. Denoten<z>≡{rj≡nzj−1 mod(d)|0<rj<d, rj+T=rj, j=1,2,3, . . . },which includes the case n=1 of cyclic (subgroup) sequence <z>. We take sets of l consecutive numbers of n<z> as T points forming a sequence{Pj(l)=(rj,rj+1, . . . ,rj+l-1)|j=1,2, . . . ,T}, in the hyper cube Cd in the 1-dimensional Euclidean space El with sides [0, d). As regards the geometry of these points, there are some crucial comprehensions that are indispensable in proceeding to full applications of spectral tests. We list their essence below.                1. The points {Pj(l)|j=1, 2, . . . } taken from n<z> in El are in a lattice (that is, occupy a portion, not necessarily all, of lattice points of a lattice) determined by d, z and l. The lattice used here agrees in definition with the one for single odd prime modulus cases in the prior art; see Detailed Description of the Invention for the definition of the lattice and its dual.        2. Spectral tests of our concern examine the configuration of lattice points in the hypercube Cd of El, not the distribution of points {Pj(l)|j=1, 2, . . . } taking their places in lattice points or seats.        3. The total number of seats in the hypercube Cd is d, irrespective of the choice of z∈Z*d or of the dimension l=2, 3, . . . .In more precise terms, we need to take all points in El that have coordinates congruent modulo d to points {Pj(l)|j=1, 2, . . . , T} noted above. This is equivalent to regard El as an l-dimensional torus with period d in all coordinates. See Nakazawa and Nakazawa (2008a,b) for plain but unabridged expositions including proofs and visual information. The total number of elements in the group Z*d is φ(d)<d. And the total number (or the period) T of elements in the coset n<z> is a divisor of φ(d) by Lagrange's theorem and T≦φ(d) is always the case. Therefore, points in El, generated by the l-tuples of <z> or n<z> in a period T leave always vacant or unoccupied lattice points in the hypercube Cd. These complications cannot be examined by spectral tests. Whatever d and z∈Z*d may be, the lattice is a definite existence, and a spectral test in El computes the largest distance λd(l)(z) between neighboring (l−1) dimensional hyperplanes that contain l−1 lattice vectors with linear independence. If z, z′ are multipliers for the same modulus d and λd(l)(z)<λd(l)(z′) holds true, the configuration of seats in El arising from z is valuated better than that given by z′. See below for more intuitive meanings of this criterion. The valuation of spectral tests has distances from the desired information on the distribution of points formed by <z> or n<z> in El. It should be stressed that the excellent performance in the spectral test is an indispensable qualification for any usable multiplier z, but not sufficient. We should be prepared for the existence of blind spots of spectral tests, and should proceed with their valuations carefully, helped by other information such as group structures of cyclic sequences or by various technical know-how that is to be accumulated in practice.        
In order to have a clear overview on spectral tests, resume the simplest case of an odd prime modulus d=p and its primitive root multiplier z. The number of lattice points is d=p, while the period T of <z> is the order of Z*p given by T=φ(p)=p−1. Points in El,{Pj(l)|j=1,2, . . . ,T}, generated by l-tuples of the sequence <z> for any dimension l=2, 3, . . . , therefore, fill the lattice points almost fully, leaving only the origin vacant. Thus, the spectral test of a primitive root multiplier of an odd prime modulus gives an almost precise valuation of the geometrical distribution of points generated by one period of <z> in El. FIGS. 2 to 4 depict visible cases of l=2 in the plane E2; the modulus is d=p=251; the multipliers are primitive roots z=162, 61 and 127, respectively; periodic sequences <z> and n<z> differ only in the choice of starting values and correspond to the same plot of T=250 points. The right of each plot shows points with normalized coordinates{(rj/d,rj+1/d)|j=1,2, . . . ,T}within a square that indicates the range −0.05≦x, y≦1.05 in E2 or xy-plane. Left plots show points with appropriately normalized coordinates{(xj/z,xj+1/z)|j=1,2, . . . ,T}corresponding to quotients that arise in the division process of the irreducible fraction n/d=1/251 to bases z=162, 61 and 127, respectively for FIGS. 2, 3 and 4. See that closeness of point configurations on the left and the right of these figures are in the order of FIG. 2, FIG. 4 and FIG. 3, in proportion to respective precision 1/z=1/162, 1/127 and 1/61. All of these figures contain the same number T=p−1=250 of generated points, or d=251 lattice points including the origin. There are wide differences in distributions of these points. If the largest distance between neighboring lattice lines (i.e. neighboring lines including the same l−1=1 linearly independent, or nonzero lattice vector in the present cases in E2) is larger, then the square can contain lesser number of those lattice lines, so that the density of points on the relevant lattice lines increases in inverse proportion. The lattice with the smallest value of its largest distance between neighboring lattice lines has the smallest density of points on lattice lines, realizing the most even or unbiased distribution of points. This will be seen clearly by comparing FIGS. 2-4. Our intuition that FIG. 4 is the worst and FIG. 2 is the best, as the distribution of points representing 2-tuples of numbers appearing uniformly and independently in sequence, is thus made quantitative in the spectral test as the largest distance λd(2)(z) of lattice lines for respective multiplier z; the smaller is λd(2)(z), the better is the multiplier z. Among all lattices in 2-dimension which are periodic with period 1 along coordinate axes and with d lattice points in the unit square, the triangular lattice realizes the smallest (hence the best) value of the largest distance between neighboring parallel lattice lines. However, in spectral tests with l-consecutive numbers generated by <z>, lattice vectors can have only rational components; see explicit forms of lattice vectors in terms of d, z and l given in Detailed Description of the Invention. Since the triangular lattice requires lattice vectors with irrational components, λd(2)(z) is invariably larger than this ideal case of the triangular lattice. The matter is the same with any other dimension l. The criterion, that the configuration is better in E2-spectral test for the lattice closer to the triangular lattice, will be utilized below as a quick and useful way of comprehension.
An epoch was marked by Fishman and Moore (1986) who initiated exhaustive spectral tests that examine all primitive roots for the prime modulus d=p=231−1. They adopted the criterion that a primitive root multiplier z is passable if λd(l)(z) is within the bound a=125% of geometrically ideal cases of the lattice (similar to the triangular lattice in 2-dimension) in all of dimensions 2≦l≦6. We note the smallest value λd(l), realized by the geometrically ideal case of lattices, for dimensions l=2, 3, . . . , 6 below citing Fishman and Moore (1986), but with a normalization that defines the lattice with period d in coordinates of El, which is adapted for the representation of lattice vectors with integer coordinates to be given later.
dimension l23456 λd(l)2−1/231/4d1/22−1/6d2/32−1/4d3/42−3/10d4/52−1/231/12d5/6Fishman and Moore (1986) note that, when the upper bound a for passable primitive root multiplier z was taken more loosely as a=133.3%, intractably many primitive roots passed the test, so that they tightened to a=125%, reducing the total number of successful primitive roots jumpingly to 414, viz. 7.74×10−7%≅1/(1.29×106) of the total number φ(231−2)=5.346×108 of primitive roots for the prime p=231−1. Later for another heavily used modulus of the type d=2i Fishman (1990) reported exhaustive spectral tests on relevant multipliers z≡5 mod(8) for d=232, and also gave tests on the modulus d=248 examining (effectively) 227 of multipliers among 245 total candidates, viz. 1/(2.6×105) of relevant generators. Remembering with dearness that Cray 2 in 1991 recorded 2 giga flops and thinking over supercomputers of today, we may contend that exhaustive spectral tests for the modulus d=248 are now realizable to obtain the period T=246. Yet this period will not suffice for simulations on supercomputers today.                Fishman and Moore (1986)/G. S. Fishman and L. R. Moore: An exhaustive analysis of multiplicative congruential random number generators with modulus 231−1. SLAM Journal on Scientific and Statistical Computing, Vol. 7 (1986), pp. 24-45.        Fishman (1990)/G. S. Fishman: Multiplicative congruential random number generators with modulus 2β: An exhaustive analysis for β=32 and a partial analysis for β=48. Mathematics of Computation, Vol. 54 (1990), pp. 331-334.        
In preceding paragraphs features of spectral tests in prior art were seen by cases in which the set of points, generated by cyclic sequences and their cosets, occupy almost all of lattice points. A significant insight on the geometry of generated points, not the geometry of lattice points, are obtained from the simplest case of composite modulus d=p1p2 formed by two distinct odd primes p1 and p2. Structures of cyclic and coset sequences <z> and n<z> for this modulus have been analyzed in Theorems A, C. Visual information will give the quickest comprehension of implications. FIG. 5A depicts consecutive 2-tuples taken from one period of the cyclic sequence <z>, emitted by a generator embodying a simple case of the present invention, obeying prescriptions of Theorem A for the case s=2. The generator has the modulus d=2867=47×61. Note that p1=47 and p2=61 giveq1=(p1−1)/2=23; q2=(p2−1)/2=30,which are coprime and realize the case (I) of Theorem C. Each square represents the range −0.05≦x, y≦1.05 again. The multiplier z isz=678≡20U1+7U2 mod(2867)where z1=20 and z2=7 are primitive roots of primes p1=47 and p2=61, respectively, both with high scores in their respective 2-dimensional spectral tests. The group Z*2867 has the order φ(2867)=φ(47)φ(61)=46×60=2760, and the period of <z> and its cosets are the largest T=2q1q2=1380=2760/2. Points on the right square represent (rj/d, rj+1/d), while those on the left square show (xj/z, xj+1/z) obtained from the numeral sequence of the irreducible fraction 1/d=1/2867=0.{dot over (x)}1x2 . . . {dot over (x)}T to base z=678. Since 1/z is sufficiently small, points on the right and the left are close to each other. FIG. 5B shows FIG. 5A with added points; the right is supplemented with points from the coset—<z>, and the left is added of points from a period of base-678 numeral sequence of the irreducible fraction −1/2867 or 2866/2867. Vacancies are seen in FIG. 5B to total to p1+p2−1=107 in number including the origin. Plots of FIG. 5B suggest a neat distribution of d=2867 lattice points in E2 on (nearly half portions of) which points of <z> and its cosets are nested. It is again tempting to use the sequence of <z> and its coset that together show the neat uniform distribution of numbers of Z*2867. But the successive appearance these sequences is the occurrence of random number sequence corresponding to 1/2867 and then the one to 2866/2867, which is highly improbable. Any concatenated sequence, n<z> and then −n<z>, should be rejected as implausible. Though the vacant seats in FIG. 5B are more conspicuous than in FIG. 2 with the sole vacancy at the origin, they are distributed uniformly so to say, and are not easy to be identified in the point sets of FIG. 5A. There is also the assurance of the shuffling by Chinese remainder theorem that the elements of the group Z*2867 are distributed in the interval [0, 2867) without accumulation. As regards the use of the sequence <z> or −<z> alone in this case, no statistical dubiousness is found in the present status of understanding.
FIG. 6 shows on the right the cyclic sequence <z> for d=2537=43×59 with z=289≡34U1+13U2 mod(2537). Here z1=34 and z2=13 are excellent primitive roots of p1=43 and p2=59, respectively. Since integersq1=(p1−1)/2=21, q2=(p2−1)/2=29are coprime, <z> has the largest period T=2q1q2=1218, and the group Z*2537 is divided into two by the cyclic subgroup <z> and its coset. However, both of q1 and q2 are odd, and the setting falls in the case of Theorem C (II) realizing −1∈<z>. Thus, only one of these halves of sequences <z> and n<z> may be used as uncorrelated random numbers. This inclusion of −1 in <z> manifests itself in the right plot of FIG. 6 as an obvious symmetry about the center M of the unit square; such a symmetry is absent in FIG. 5A. In fact, the following general inference holds true:                Lemma D Let there be given a modulus d>4 which may be composite. Take any z, n∈Z*d and an arbitrary integer l≧1. In the l-dimensional Euclidean space El denote Sl for the set of points with coordinates congruent modulo d to consecutive l-tuples of integers in the coset sequence n<z> in Z*d,Sl:={Pj(l)≡(nzj,nzj+1, . . . ,nzj+l-1)mod(d) j=0,1,2, . . . }.Let S′l and Sl be point sets in El symmetric to Sl with respect to the point Md=(d/2, d/2, . . . , d/2) and to the origin O, respectively.        (I) If −1∈<z> is true, there holds Sl=S′l. Namely, the point set formed by l consecutive numbers of n<z> defined modulo d is symmetric about the point Md, irrespective of the choice of n∈Z*d including the case n=1.        (II) If −1 ∉<z> is the case, then Sl∩S′l=φ holds true. Namely, the point set formed by l consecutive numbers of n<z> defined modulo d is asymmetrical about the point Md, irrespective of the choice of n∈Z*d including the case n=1.(Proof) We prove the case l=1; this will give indications of the proof for l≧2 given in Nakazawa and Nakazawa (2008a,b), where the case l=1 was not discussed. Integers of n<z> form a point set S on the line E1, to be called x-axis in this proof. Since n<z> is defined modulo d the set S is invariant under d-translations. Let S′ and S denote point sets on the x-axis which are symmetric to S about points x=d/2 and x=0, respectively; these are also d-translation invariant, and S may be denoted S=−S. The following is obvious:S′={−x+d|x∈S}={−x|x∈S}= S. If −1∈<z> holds true, then—<z> is a sequential shift of <z>, so that S=S holds as sets. Thus we obtain S′= S=S. If −1 ∉<z> is the case, then <z> is a proper subgroup of Z*d, and disjoint as a set with the coset—<z>, implying S′∩S= S∩S=φ. ▪FIGS. 2-4 show primitive root multipliers for odd prime moduluses, and cyclic sequences <z> therein all include −1. Accordingly, right plots show point sets that are symmetric w.r.t. the center of the square. Similarly, FIG. 6 shows the composite modulus case with −1 included in <z>. The right plot is again symmetric w.r.t. the center of the square. In contrast, FIG. 5A stands for the case that <z> does not include −1, and its right plot shows no conceivable symmetry. Cases (I) and (II) in Lemma D give the dichotomy, so that we may state:        Corollary to Lemma D Let z be a multiplier taken in the group Z*d of an arbitrary modulus d>4, and let Sl be the set of points in El formed by the consecutive l-tuples of a coset sequence n<z> for an arbitrary element n in Z*d. Inclusion of −1 in <z> is equivalent to the symmetry of Sl w.r.t. the center Md of the hyper cube of sides [0, d) for any one of dimensions l=1, 2, . . . .Inclusion of −1 in <z>, or the symmetry of point sets formed by l-tuples of integers of n<z> in El, is in fact a more serious symptom than apparent symmetry of formed point sets; the second half of generated points take symmetric positions reproducing the order of appearance of their correspondents in the first half.        
The present invention chooses settings of cyclic sequence <z> in the reduced residue class group Z*d for a composite modulus d, where −1 is excluded from <z> by design. An example was shown in FIG. 5A. As another we show in FIG. 1 the result for the settingd=2867=47×61, z=1298≡29U1+17U2 mod(2867),where z1=29 and z2=17 are excellent primitive roots for p1=47 and p2=61. The right plot of two-dimensional spectral test shows nearly triangular lattice. Comparisons with FIG. 5A will impress us that the multiplier z=678 for FIG. 5A will also be excellent in the two-dimensional spectral test. This is in fact the case. However, there can be a large difference in higher dimensions. The following table gives results of their spectral tests for l=2, 3 as comparisons to λd(l) given by geometrically ideal cases, indicating how can differences be for choices of the multiplier z.
2-nd degree3-rd degreemultiplier zscorescorez = 1298 ≡ 29U1 + 17U2 mod (2867)104.96%138.27%z = 678 ≡ 20U1 + 7U2 mod (2867)105.67%196.36%
We saw what to be seen. We make use of a group Z*d with a composite modulus d in order to exploit the noted improvements given by the shuffling, and also to exclude −1 from the cyclic sequence <z>. The group Z*d for any modulus d has integers in [0, d) as its elements, with vacancies that may also be said to be distributed evenly therein. We choose parameters so as for these vacancies to be small in number, and sequences <z> and n<z> will sweep as large a portion of Z*d as possible, for uniformity and for the strength of spectral tests. Lagrange's theorem restrains that the (sub)group <z> may occupy only an integral portion of Z*d. And a composite modulus group cannot be cyclic. The largest portion that <z> can occupy in Z*d is thus the half, and the structure of the period shown in Theorem A rules that the modulus should then comprise s=2 odd distinct primes. The factor d0 of 1 or 2 in the modulus has no effect on the period and its participation in the shuffling is merely as a constant factor for d0=2:nzj≡d+(d+1)(n1z1j+n2z2j+ . . . +nszsj)mod(2d),so that we discard d0 altogether. We are thus left with the form d=d1d2, with d1=p1i1 and d2=p2i2 for odd primes p1 and p2.
Component groups Z*pkik for k=1, 2 have generating elements with periodsTk=φ(pkik)=pkik(1−1/pk)≅pkik, k=1,2.For the shuffling in <z> and its coset to be effective, we require T1≅T2≅t. Then the period of <z> is T=T1T2≅t2. Now the computability of spectral tests puts a practical limit on T or on t≅T1/2. Component moduluses dk=pkik (k=1, 2) should thus be determined by the requirementt≅T1/2≅pkik, pk≅t1/ik, k=1,2.Therefore, the ratio V of the number of voids, i.e. integers in [0, d) left unoccupied by the elements of Z*d, to the total number d of lattice points isV={d−φ(d)}/d=1−(1−1/p1)(1−1/p2)≅1/p1+1/p2≅t−1/i1+t−1/i2.This factor for a fixed t is minimized at i1=i2=1. The best setting for our aim is thus provided by moduluses of the form d=p1p2 comprising two distinct odd primes of similar magnitudes, on the premise that we exploit the largest period setting of Theorem A.
Summarizingly, we have first to determine the computable largest period T=t2. Then we choose distinct odd primes, p1, p2≅t; the choice should obey prescriptions of Theorem A and Theorem C that integers q1=(p1−1)/2 and q2=(p2−1)/2 are coprime and have different parity, one even and the other odd. The multiplier z is then formed by primitive roots z1, z2 of primes p1, p2, respectively, by congruence relations z≡z1 mod(p1) and z≡z2 mod(p2).